Suppose that I want to simulate a gas, i.e. a collection of particles in a two-dimensional box that move around. A particle has a position, which is a 2D vector, a velocity, which is another 2D vector, and a mass, that is a real number. How can I represent this in Julia?
As we have started to see, the basic types in Julia, such as Complex, are defined in Julia code.
They are treated on the same footing as types defined by the user.
Composite types in Julia are just boxes that hold data and provide constructor functions for making new objects of that type. They do not contain methods (functions) (as would be the case in many other object-oriented languages); rather, the methods are defined using multiple dispatch.
Currently, there are no fixed-size vectors in base Julia (although these are available, for example, in the ImmutableArrays and FixedSizeVectors packages). For many applications, fixed-size vectors are useful, for example for representing positions of particles moving in 2D or 3D space. We can define our own type that may prove to be more efficient than Julia's Array type.
The basic syntax for defining a composite type is
In [1]:
type Vector2D
x::Float64
y::Float64
end
The double colon (::) is a type annotation that specifies the type of x and y.
Let's check which methods have been defined:
In [2]:
methods(Vector2D)
Out[2]:
In [3]:
Vector2D(3.5, 4.5)
Out[3]:
In [4]:
Vector2D(3, 4)
Out[4]:
In [5]:
v = Vector2D(3, 4)
w = Vector2D(5, 6)
Out[5]:
What happens if we try to add two Vector2Ds?
In [7]:
v + w
Julia confirms that it has no idea how to do so. So let's just try defining it:
In [8]:
+(a::Vector2D, b::Vector2D) = Vector2D(a.x+b.x, a.y+b.y)
Out[8]:
In [9]:
v + w
Out[9]:
Suddenly, it works! This is Julia's version of operator overloading -- just add a new method to the relevant operator that acts on objects of your new type.
In [ ]:
dot(a::Vector2D, b::Vector2D) = a.x*b.x + a.y+b.y
Julia has certain Unicode operators defined, for example \cdot<TAB>:
In [ ]:
⋅
We see that this operator is just an alias for the dot function, but that is parsed as in infix operator. So we can now write
In [ ]:
v ⋅ w
It doesn't work. What's the matter? The signal was given by Julia when we defined dot: it said there was only 1 method. But \cdot<TAB> has 7 methods. So they are different functions. The solution is that the name dot refers to Julia's built-in dot function, which is defined in Julia's standard library, called Base. So the true name of the function is Base.dot:
In [ ]:
Base.dot
We thus need to extend this function instead:
In [ ]:
Base.dot(a::Vector2D, b::Vector2D) = a.x*b.x + a.y+b.y
In [ ]:
v ⋅ w
Julia provides an import keyword to make this simpler when we are defining different methods for the function:
In [ ]:
import Base.dot
This must be done before defining dot. To clear our current namespace we can use workspace:
In [ ]:
workspace()
In [ ]:
dot
In [ ]:
type Vector2D
x::Float64
y::Float64
end
In [ ]:
import Base.dot
In [ ]:
dot(a::Vector2D, b::Vector2D) = a.x*b.x + a.y+b.y
In [ ]:
v = Vector2D(1, 2)
w = Vector2D(3, 4)
In [ ]:
v ⋅ w
We may wish to define different kinds of Vector2D that contain different types. We can define them all at once by using a type parameter:
In [ ]:
workspace()
In [ ]:
immutable Vector2D{T}
x::T
y::T
end
We have replaced the keyword type by immutable; this permits the compiler to store the object more efficiently. The {T} indicates that T can now be any type, and x and y are both of the same type T.
Exercise: Create vectors of different types
Functions may also be parametrised:
In [ ]:
+{T}(a::Vector2D{T}, b::Vector2D{T}) = Vector2D{T}(a.x+b.x, a.y+b.y)
Exercise: Define a particle type and create a couple of particles. Define a function move! that moves a particle in its current direction for a small time $\delta t$.
Exercise: Define a gas that contains a certain number of particles, $N$, and the particles themselves. Define move! for the gas. (A gas should also have a box in which it lives -- exercise.)
As in other languages, a type's behaviour is determined both by its fields and by the methods that act on the types; the only difference, once again, is that in Julia the methods are defined outside the definition of the type itself.
We will now define another type that has two fields but that behaves very differently, to permit automatic (or algorithmic) differentiation.
The idea is that to differentiation a complicated function, say
$f(x) := \sin[(x-1)(x^2+2)],$
the function is decomposed into elementary functions that we know how to differentiate.
For example, suppose that we know how to differentiate the functions $g(x) := x-1$ and $h(x) := x^2 + 2$ at $x=a$. Then the derivative of their product is given by
$[g \cdot h]'(a) = g(a) h'(a) + g'(a) h(a)$.
Similarly, to differentiate any function at $a$, it is sufficient to know the value and the derivative of each subpart of the function at the same point $a$.
We thus represent a function $f$ by the pair $(f(a), f'(a))$, which we can regard as a mathematical object called a jet), i.e. the collection [equivalence class] of all functions that have the same value and the same derivative at the point $a$.
In [ ]:
immutable Jet
value::Float64
deriv::Float64
end
The operations on Jets should correspond to the operations on those functions to give new functions. For example, we have
In [ ]:
+(f::Jet, g::Jet) = Jet(f.value+g.value, f.deriv+g.deriv)
*(f::Jet, g::Jet) = Jet(f.value*g.value, f.value*g.deriv + f.deriv*g.value)
A constant function $c(x) := c$ for all $x$ has derivative $0$, so c corresponds to the jet Jet(c, 0).
The identity function $\mathrm{id}(x)$ has the value $id(a)=a$ at the point $x=a$ and derivative $1$, and so corresponds to the Jet Jet(a, 1).
We can now do simple derivatives. Let's calculate
In [ ]:
a = 3
xx = Jet(a, 1)
f(x) = x*x + x # we haven't yet defined ^ for a Jet!
f′(x) = 2x + 1
In [ ]:
@show f(a)
@show f′(a)
@show f(xx)
This correctly calculates, automatically, the value and derivative of the function at $a$!
Suppose we now try to do
In [ ]:
xx + 3
Julia does not know how to add a jet and an integer -- because we defined addition only for two jets.
One (painful) solution would be to explicitly define +(a::Jet, b::Number), +(a::Number, b::Jet), and the same for all the other operations.
Julia provides a better solution: a system for promotion and conversion of types from one type to another. Nothing is "automatic", though, it's all coded in rules. For example, what does Julia do when we try to add an integer and a float?
In [ ]:
@which 3 + 4.5
We see that it calls a very generic function that adds two Numbers. We copy here the code for convenience:
In [ ]:
+(x::Number, y::Number) = +(promote(x,y)...)
What does this call to promote do?
In [ ]:
promote(3, 4.5)
It returns both numbers, converted to the same type. How does it work?
In [ ]:
@which promote(3, 4.5)
That code is
In [ ]:
function promote{T,S}(x::T, y::S)
(convert(promote_type(T,S),x), convert(promote_type(T,S),y))
end
It uses the function convert that converts between types, e.g.
In [ ]:
convert(Float64, 3)
In [ ]:
@which convert(Float64, 3)
That code calls out the sitofp LLVM intrinsic function that converts a signed integer to a floating point number.
What does promote_type do?
In [ ]:
promote_type(Float64, Int)
In [ ]:
@which promote_type(Float64, Int)
In [ ]:
promote_type{T,S}(::Type{T}, ::Type{S}) =
promote_result(T, S, promote_rule(T,S), promote_rule(S,T))
And finally,
In [ ]:
promote_rule(Float64, Int64)
In [ ]:
@which promote_rule(Float64, Int64)
This specifies what the result of the promotion of those two types should be.
We can "hook in" to this mechanism if we think of a Jet as a type of Number. We can specify that Jet is a new type in the type hierarchy below Number as follows:
In [ ]:
workspace()
In [ ]:
immutable Jet <: Number
value::Float64
deriv::Float64
end
In [ ]:
Base.convert(::Type{Jet}, c::Real) = Jet(c, 0)
Base.promote_rule{S<:Real}(::Type{Jet}, ::Type{S}) = Jet
In [ ]:
+(f::Jet, g::Jet) = Jet(f.value+g.value, f.deriv+g.deriv)
*(f::Jet, g::Jet) = Jet(f.value*g.value, f.value*g.deriv + f.deriv*g.value)
Now operations with constants work "magically" (without any magic at all, i.e. with everything being completely explicit and traceable in the code):
In [ ]:
j = Jet(3, 1)
In [ ]:
j + 2
In [ ]:
j * 2
In [ ]:
2 * j
So now we can differentiate more complicated functions:
In [ ]:
f(x) = 3x^2 + 2.3x
In [ ]:
a = 3
xx = Jet(a, 1)
f(xx)
Exercise: Define ^ for jets and check that you can differentiate polynomials. Define exp for a jet and differentiate functions involving exponentials.
Note that several implementations of automatic differentiation are available.
An implementation along the lines presented here is available in ValidatedNumerics.jl.